!**************************************
MODULE DATA_ARRAY
    real,dimension(1999,1999) :: Rij_n,Rij
    integer :: pairn,pn
    integer,dimension(1999) :: pi,pj
END MODULE DATA_ARRAY
!**************************************


program main
use DATA_ARRAY
implicit none

 character(len=99)buffer,maplist1,trjpdb,output
 character(len=99)buffer1,maplist2
 character(len=20),dimension(20) :: col
integer i,j,k,ios,totframe,nargs,connum1,connum2,resi,resj,l,m,n,transNum
integer,parameter :: Contactnumber=1999
integer,parameter :: closeca=196
integer :: c1(Contactnumber),c2(Contactnumber)
integer :: c3(Contactnumber),c4(Contactnumber)
real :: Q(Contactnumber),Qnative,dist,P(closeca,closeca,2),Qf,Qplu,d(71,71)
real,parameter :: cutoff1=6.0
real,parameter :: cutoff2=8.5
real :: x(closeca)
real :: y(closeca)
real :: z(closeca)
!real :: d1(Contactnumber),d2(Contactnumber)
!real :: R0(Contactnumber),R1(Contactnumber)
real,parameter :: MULTI=1.2

 nargs=iargc()
   if(nargs==4)then
      call getarg(1,maplist1)
      call getarg(2,maplist2)
      call getarg(3,trjpdb)
      call getarg(4,output)
   else
      write(*,*) "Wrong Input!"
      write(*,*) "#Usage: /xx.x  contactmaplist  trjpdb   outputfile"
      write(*,*) "Usage: /xx.x  binding73.list  folding_csu.list totalpdb.pdb   distance.list"
      write(*,*) "Usage: /xx.x  inter-pairs.list folding.list traj_tmp.pdb   distance.list"
      goto 1000
   endif

ios=0
i=0
open(10,file=maplist1,status="old")
do while(.true.)
 read(10,"(A)",iostat=ios)buffer
 if(ios/=0)exit
 i=i+1
 read(buffer,*)(col(k),k=1,3)
 read(col(1),"(I5)") resi
 read(col(2),"(I5)") resj
 read(col(3),"(F10.7)") dist
 c1(i)=resi
 c2(i)=resj
 Rij_n(resi,resj)=dist
end do
close(10)

 connum1=i
!    write(*,*) connum1
ios=0
i=0
open(11,file=maplist2,status="old")
do while(.true.)
 read(11,"(A)",iostat=ios)buffer
 if(ios/=0)exit
 i=i+1
 read(buffer,*)(col(k),k=1,3)
 read(col(1),"(I5)") resi
 read(col(2),"(I5)") resj
 read(col(3),"(F10.7)") dist
 c3(i)=resi
 c4(i)=resj
 Rij_n(resi,resj)=dist
end do
close(11)

 connum2=i
!    write(*,*) connum2

ios=0
k=0
transNum=0
Qnative=0
Qf=0
open(20,file=trjpdb,status="old")
open(30,file=output,status="replace")
!open(35,file="QiQf.dat",status="replace")
open(40,file="transition.list",status="replace")
!open(50,file="nonnative.list",status="replace")
do while(.true.)
 read(20,"(A)",iostat=ios)buffer
  if(ios/=0)exit
  if(buffer(1:5)=="MODEL")then
   i=0
   k=k+1
!    if(mod(k,1000)==0)then
!      write(*,*) k
!    endif
!   Qnative=0
!   Qf=0
    do while(.true.)
     read(20,"(A)")buffer
      if(buffer(1:3)=="TER")exit
      if(buffer(14:14)=="C")then
        i=i+1
       read(buffer(32:38),"(F)") x(i)
       read(buffer(40:46),"(F)") y(i)
       read(buffer(48:54),"(F)") z(i)
      end if
    end do
    do j=1,connum1
     Rij(c1(j),c2(j))=sqrt((x(c2(j))-x(c1(j)))**2+(y(c2(j))-y(c1(j)))**2+(z(c2(j))-z(c1(j)))**2)
!      if(d1(j)<MULTI*R0(j))then
!          Qnative=Qnative+1.0
!      endif
    end do
    do j=1,connum2
     Rij(c3(j),c4(j))=sqrt((x(c4(j))-x(c3(j)))**2+(y(c4(j))-y(c3(j)))**2+(z(c4(j))-z(c3(j)))**2)
    end do
!===============================================================
!  calculate RC:Qb,Qf,and determine the states according to RCs
!===============================================================
    Qplu=0.0
    pi=c1
    pj=c2
    pn=connum1
    call qcal(Qplu)
    Qnative=Qplu/pn

    Qplu=0.0
    pi=c3
    pj=c4
    pn=connum2
    call qcal(Qplu)
    Qf=Qplu/pn

!       write(35,"(I,3X,2F10.7)") k,Qnative,Qf
!       if(Qnative<0.9 .and. Qnative>0.5 .and. Qf<0.5 .and. Qf>0.2)then
!       if(Qnative<=0.5 .and. Qf<0.5)then
       if(Qnative==0.0 )then
!       if((Qnative>0.0 .and. Qf<0.5) .or. (Qnative<=0.7 .and. Qf>=0.5))then
!    do j=1,connum1
!!     d1(j)=sqrt((x(c2(j))-x(c1(j)))**2+(y(c2(j))-y(c1(j)))**2+(z(c2(j))-z(c1(j)))**2)
!      if(Rij(c1(j),c2(j))<MULTI*Rij_n(c1(j),c2(j)))then
!          Q(j)=Q(j)+1
!      endif
!    end do
       transNum=transNum+1
!	do m=1,71
!		do n=72,closeca
!		P(m,n,1)=sqrt((x(m)-x(n))**2+(y(m)-y(n))**2+(z(m)-z(n))**2)
!		if(P(m,n,1)<=cutoff1)then
!		P(m,n,2)=P(m,n,2)+1.0
!		elseif(P(m,n,1)<=cutoff2 .and. P(m,n,1)>cutoff1)then
!		P(m,n,2)=P(m,n,2)+0.5
!		endif
!		enddo
!	enddo
	do m=1,71
		do n=1,71
		 d(m,n)=d(m,n)+sqrt((x(m)-x(n))**2+(y(m)-y(n))**2+(z(m)-z(n))**2)
		enddo
	enddo
       write(40,"(I,3X,F15.9,3X,F15.9)") k,Qnative,Qf
       endif
   end if
enddo
  totframe=k
!   write(*,*) "Total frames :  ",totframe
!   write(*,*) "Total transition frames :  ",transNum
!!       write(30,"(I,3X,F15.9)") k,Qnative/(connum*totframe)
!! !      write(30,*) "##Atomi  Atomj  Qij  Resi  Resj##"
!    do j=1,connum1
!       write(30,"(I3,2X,I3,3X,F7.3)") c1(j),c2(j),Q(j)/transNum
!    enddo
!	do m=1,71
!		do n=72,closeca
!		if(P(m,n,2)/=0)then
!		write(50,"(I5,3X,I5,3X,F7.3)") m,n,P(m,n,2)/transNum
!		endif
!		enddo
!	enddo
	do m=1,71
		do n=1,71
		 write(30,*) m,n,d(m,n)/transNum
		enddo
	enddo
close(20)
close(30)
!close(35)
close(40)
!close(50)
1000 end program main

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine qcal(Qplu)
use DATA_ARRAY
implicit none
real,intent(out) :: Qplu
integer :: ipair,i,j,npair
real :: R0,R1
real,parameter :: RATE=1.2

Qplu=0.0
npair=pn

 do ipair=1,npair
   i=pi(ipair)
   j=pj(ipair)
   R0=Rij_n(i,j)
   R1=Rij(i,j)
	if( R1 < (RATE * R0) )then
	 Qplu=Qplu+1
	endif
 enddo
end subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!subroutine qcal(Qplu)
!use DATA_ARRAY
!implicit none
!real,intent(out) :: Qplu
!integer :: ipair,i,j,npair
!real :: R0,R1,D0
!real,parameter :: m=20,n=10,l=0.30
!
!Qplu=0.0
!npair=pn
!
! do ipair=1,npair
!   i=pi(ipair)
!   j=pj(ipair)
!   R0=Rij_n(i,j)
!   R1=Rij(i,j)
!   D0=l*R0
!   Qplu = Qplu + (1-((R1-D0)/R0)**n) / (1-((R1-D0)/R0)**m)
! enddo
!end subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
